Elliptic Grid Generation with B-Spline Collocation

We revisit the lassi al te hnique of ellipti grid generation with harmoni mappings. For the determination of the ontrol fun tions we use the framework developed by Spekreijse [1℄. However, instead with nite di eren es we dis retize the underlying partial di erential equation with a B-spline ollo ation method in order to work dire tly with the native data representations of our CAGD system. This way we an make use of the sparsity and a ura y of the B-spline boundary representations and guarantee the geometri onsisten y of our CAD models. In this paper we will summarize the underlying algorithms and present some rst appli ation examples. Introduction In the ontext of the development of a new adaptive Navier-Stokes solver quadflow whi h aims at the simulation of uid-stru ture intera tion at airplane wings, f. [3℄, a new grid generation module has been implemented whi h is based on the representation of the geometry with parametri mappings. As in many ommon CAD systems free form urves and surfa es are represented by B-splines. From this parametri representations adaptive grids are omputed by fun tion evaluation. Con retely, within this system urves are represented by B-splines of the form x(u) = N Xi=0 piNi;p;U (u); (1) planar grids and surfa es are modeled by bivariate B-spline tensor produ ts x(u; v) = N Xi=0 M Xj=0 pi;j Ni;p;U (u)Nj;q;V (v) (2) and volume grids are represented by trivariate mappings x(u; v; w) = N;M;L X i;j;k=0pi;j;kNi;p;U (u)Nj;q;V (v)Nk;r;W (w): (3) Here U; V;W are non-de reasing and non-stationary knot sequen es, i.e., U = (ui)N+p 2 i=0 : u0 u1 : : : uN+p 2; ui+p < ui; (4) and Ni;p;U denotes the i-th normalized B-spline-fun tion of order p de ned by the re ursion Ni;1;U (u) = (1 if ui u < ui+1 0 otherwise ; (5) Ni;p;U (u) = u ui ui+p 1 uiNi;p 1(u) + ui+p u ui+p ui+1Ni+1;p 1(u): (6) B-spline urves are pie ewise polynomials of degree p 1. Usually we hoose p = 4, i.e. ubi splines, and knot sequen es with p-fold knots at the interval ends. This has the advantage, that the knot sequen e interval oin ides with the parameter interval of the urve, that the rst and last ontrol point oin ide with the start and end point of the urve, and that the rst and last span of the ontrol polygon are tangential to the urve at the start and end point of the urve. In many pra ti al ases non-uniform knot sequen es are onstru ted, for instan e in the ourse of an adaptive approximation or interpolation pro edure. Multiple interior knots an be used to model non-smooth features of an obje t. Grid Generation Equations In order to generate smooth grids, all the ommon te hniques of stru tured grid generation, namely trans nite interpolation and methods based on the solution of partial di erential equations are applied. In this paper we want to integrate an ellipti grid generation te hnique into our CAGD system. Our hoi e was for Spekreijse's approa h whi h an be very brie y summarized as follows: let x(s) be a harmoni mapping from the d-dimensional parameter spa e P onto the physi al domain C and s( ) be a soalled ontrol mapping from the omputational domain C onto the parameter domain P . Then the omposite mapping x(s( )) : C ! D (7) ful lls a di erential equation of the form L(x) = d X i;j=1 gij 2x i j + d X k=1Pk x k = 0 (8) where Pk = d X i;j=1 J2gijP k ij ; (9) J = det x0( ) is the Ja obian of the omposite mapping, the gij and gij are the ovariant and ontravariant metri tensors de ned by gij = x i x j ; d Xk=1 gikgkj = Æij ; (10) the P k ij are the omponents of the ve tor Pij = T 1 2s i j ; (11) and T = s0( ) is the Ja obian matrix of the ontrol mapping. The aim of this paper is to solve this PDE with a B-spline ollo ation method. B-Spline Collocation The general idea of ollo ation is to determine a fun tion, so that it exa tly satis es the di erential equation at ertain points, the ollo ation points. In a way ollo ation is similar to interpolation, but in ontrast to interpolation we do not mat h fun tion values but ertain ombination of fun tion and derivative values. In order to simplify the notation we on entrate onto the bivariate ase from now on and denote the Cartesian oordinates of the omputational domain, the unit square, with = (u; v) and of the parameter domain with s = (s; t). Hen e, we sear h a fun tion of the form (2) whi h ful lls Lx(ûi; v̂j) = 0; i = 1; : : : ; N 1; j = 1; : : : ;M 1 (12) at ertain ollo ation points ûi, v̂j and Diri hlet boundary onditions for the ontrol points p0;j , pN;j , j = 0; : : : ;M and pi;0, pi;M , i = 0; : : : ; N . The task is now to hoose appropriate ollo ation points for the on guration at hand. The most popular B-spline ollo ation s heme is Gauÿ ollo ation with ubi Hermite-splines. Its appli ation to ellipti grid generation has already been investigated by Manke [2℄. Here in ea h knot interval the ollo ation points are the abs issae of Gaussian quadrature rules. An obvious aveat is that in a prepro essing step one has to make all interior knots two-fold by knot-insertion, and that therefore the resulting grid will be only C1. Our preferen e is for a s heme that works for splines with arbitrary knot sequen es and uses the Greville abs issae, whi h are de ned by ui = i+p X k=i+1 uk; (13) as ollo ation points. This hoi e is motivated by the S hoenberg-Whitney theorem, see referen e [4℄, whi h says that the interpolation problem x(ûi) = fi is well posed if, and only if, every ûi lies in the support of the i th B-spline fun tion, i.e., if Ni(ûi) > 0. As one an easily verify, the Greville abs issae always give a set of as many distin t points, as the spline has ontrol points and they ful ll the onditions of the S hoenberg-Whitney theorem. This ollo ation s heme frees us from the ne essity to insert additional knots in our tensor produ t representations. A disadvantage, however, is that ollo ation at the Greville abs issae does not have the optimal onsisten y order. The S hoenberg-Whitney theorem is also the reason why we do not use a standard nite di eren e ode followed by an interpolation algorithm in order to generate ellipti B-spline grids: a typi al nite di eren e ode is based on the assumption that the grid points xi;j are numeri al approximations of regular spa ed values x(ihu; jhv). However, depending on the stru ture of the underlying spline it ould be ome ne essary to work with unevenly spa ed grids in order to ful ll the stipulations of the S hoenbergWhitney theorem during the interpolation pro ess. Application Example The afore-mentioned ollo ation s hemes have been implemented and tested for planar grids, surfa es and volume grids. In order to solve the systems we just follow the standard approa h and use a xed point iteration, freezing the metri oe ients in Equation (8) in order to get a linear system in every single iteration. Then we apply the ollo ation s heme to the linearized equations. The arising sparse linear systems are solved with a dire t solver. This kind of implementation is not well suited to solve big systems with maximum e ien y, but the aim of the urrent study was not to ompare the e ien y of the implementation (this has already be done in [2℄), but to study the prin ipal validity of the method. As a rst appli ation we present the grid in a blo k that is taken from a grid for a dual-bell on guration, see Figure 1. The boundaries are approximatively parameterized by ar length, so that we an use the identi mapping as ontrol mapping. Hen e, all ontrol fun tions P k ij are zero and the resulting grid mapping is harmoni . However, the spa ing of the ontrol points, that an be seen in the upper plot is rather irregular. This irregularity stems from an adaptive B-spline approximation algorithm whi h tried to resolve the di erent features of the nozzle ontour and from the ne essity to mutually insert the knots whi h are not present in the representation of the opposite boundary in order to build a tensor produ t. However, the resulting numeri al grid, whi h is omputed by evaluation of the B-spline fun tion has the desired smoothness properties. Figure 1. Control points and harmoni mesh for the dual bell. Boundary Orthogonality In Spekreijse's approa h there remains the problem to determine suitable ontrol fun tions in order to in orporate desired features into the grids. In order to nd a ontrol mapping that ensures boundary orthogonality Spekreijse proposes to pro eed as follows. Let us assume that a folding-free grid x( ) is already available. This grid may be, for instan e, a trans nite interpolant or the solution of the purely harmoni grid generation system. Then in the rst step we solve the transformed Lapla e equation div(A grad s) = 0 (14) where A = J g11 g12 g12 g22 = 1 J g22 g12 g12 g11 : (15) This equation is supplied by mixed Diri hlet and Neumann boundary onditions, in parti ular we require s= n = 0 at the boundaries x(u; 0) and x(u; 1) and t= n = 0 at the boundaries x(0; v) and x(1; v) of the physial domain. The solution of this problem gives us a one-to-one boundary mapping C ! P . In the se ond step we omplete this boundary mapping to a suitable ontrol mapping that ful lls the orthogonality onditions t= u = 0 along the boundaries u = 0 and u = 1 and s= v = 0 along the boundaries v = 0 and v = 1 using an algebrai grid generation method. For the details of this method we have to refer to [1℄. Whereas the dis retization of Equation 8 by ollo ation is straightforward it is more onvenient to dis retize equation 14 with a nite volume method. For this we observe that for any ontrol volume in the omputation domain the equation Z (grads; An) d = 0 (16) holds. Of ourse, as ontrol volumes we will hoose intervals of the form [ui; ui+1℄ [vj ; vj+1℄. Again we want to represent the ontrol mapping as tensor produ t B-spline. Therefore, the integral over the boundary of the ontrol volume is omposed of segments of the form Z ui+1 ui (grads; An) d = Z ui+1 ui su sv ; A 01 d =Xi;j pij Nj(v) Z ui+1 ui N 0 i(u)g12 J du+N 0 j(v) Z u2 u1 Ni(u)g11 J du and Z vj+1 vj (grads; An) d = Z vj+1 vj su sv ; A 10 d =Xi;j pij Ni(u) Z vj+1 vj N 0 j(v)g12 J dv +N 0 i(u) Z vj+1 vj Nj(v)g22 J du : The integrals in the bra kets enter the matrix of the dis retized problem and an heaply be evaluated by quadrature formulas. In order to get as many equations as ontrol points we enter the ontrol volumina around the Greville abs issae by hoosing ui = 1 2( ui 1 + ui); i = 1; : : : ; N 1; u0 = 0; uN = 1 vi = 1 2( vj 1 + vj); i = 1; : : : ; N 1; v0 = 0; vM = 1 (17) The boundary ondition s= n transforms to (grads; An) = 0 at the orresponding boundary in C, so that the dis retization at boundary grid points is also straightforward. Figure 2 shows the ontrol points and the smooth evaluation of the resulting orthogonal grid, Figure 3 shows the orresponding ontrol mapping and a detail view at the nozzle throat. Figure 2. Control points and orthogonal mesh for the dual bell. Figure 3. Control Map and Detail View Convergence and Stability Matters The above example shows that in prin iple the ollo ation method presents a viable method to integrate ellipti grid generation strategy into a spline based CAGD system. However, it turns out that there are also some ompli ations. First signs for these problems already reveal themselves in the following onvergen e study. Figure 4. Cosine Test ase We onsider a simple re tangle where the lower side has been repla ed by a osine-like ar , see Figure 4, and ompare the L2-residual of the fully onverged solutions, i.e., the solution we get when the xed point iteration does not improve the solution any more. Gauÿ Greville N rN := kL(x)k2 rN=2=rN rN := kL(x)k2 rN=2=rN 10 4.832e+2 4.757e+2 20 3.138e+2 1.54 2.916e+2 1.63 40 1.233e+1 2.55 1.480e+1 1.97 80 6.146 2.01 5.939 2.49 160 2.570 2.39 1.988 2.99 Figure 5. Convergen e behavior of the osine test ase. First of all we observe that the Gauÿ ollo ation s heme does not produ e better onvergen e faster than the Greville ollo ation s heme, although theory predi ts fourth order onvergen e for the former and only se ond order for the latter. (We an indeed observe these rates for the linear Lapla e equation with our implementation). The Hermite s heme even onverges slower than the Greville s heme, what somehow justi es our preferen e forthe Greville s heme. One reason for this might be the very bad ondi-tion of the ollo ation matri es whi h typi ally goes with No where o isthe order of the s heme. However, from the geometri point of view theresult is not entirely satisfa tory either. As is well known, harmoni gridgeneration systems have the tenden y to push away the grid lines fromon ave boundaries. Espe ially if one tries to apply the ollo ation s hemewith very oarse ontrol nets this tenden y seems to be even a entuated.Figure 4, for example, shows the resulting 10 10 ontrol point grid andthe orresponding grid evaluation. This defe t disappears in the ourse ofextensive grid re nement, but only very slowly. For instan e, in the evenmore extreme example of Figure 6 one needs more than 160 ontrol pointsin ea h oordinate dire tion before the grid lines start to onverge towardsthe boundary. The Figure itself shows the result of the dis retization basedon 40 40 ontrol points.Figure 6. Dis retization with 40 40 ontrol pointsAt this pla e it is interesting to note, that the solution of the transformedLapla e equation (14) an also be used to ompute for any given grid map-ping x(u; v) a orresponding ontrol mapping s(u; v). This an be done byrepla ing the Neumann boundary onditions by standard Diri hlet ondi-tions. Figure 7 shows the ontrol mapping that orresponds to the grids inFigure 6. The aberration of this grid from the identi mapping is obviouslyaused by the dis retization error.ConclusionThe here proposed ollo ation s heme presents a useful method to real-ize ellipti grid generation methods dire tly on B-spline representations.However, one needs additional measures to ontrol the boundary spa ing. Figure 7. Control mapping orresponding to the grid in Figure 6.Following Spekreijse this would a ord the omputation of appropriate on-trol mappings via solution of the biharmoni equation. This, however, hasnot yet been implemented for B-spline grids. Sin e in [2℄ Manke does notreport similar problems, it might well be that methods whi h ompute theontrol fun tions iteratively are a reasonable alternative in this ontext,too.AcknowledgmentsThis work has been and performed with funding by the Deuts he Fors hungs-gemeins haft in the Collaborative Resear h Center SFB401 "Flow Modu-lation and Fluid-Stru ture Intera tion at Airplane Wings" of the RWTHAa hen, University of Te hnology, Aa hen, Germany.References[1℄ Stefan Spekreijse. Ellipti Grid Generation Based on Lapla e-Equationsand Algebrai Transformations. J. Comp. Phys., 118:38 61, 1995.[2℄ J. W. Manke. A Tensor Produ t B-Spline Method for Numeri al GridGeneration. J. Comp. Phys., 108:15 26, 1993.[3℄ Frank Bramkamp, Philipp Lamby, and Siegfried Müller. An adaptivemultis ale nite volume solver for unsteady and steady state ow om-putations. J. Comp. Phys., 197:460 490, 2004.[4℄ Carl de Boor. A Pra ti al Guide to Splines. Springer, 1978.